Bray-curtis
Proximal small intestine
set.seed(123)
# subset dateframe
df_abundance <- subset(df_ST1T, select = -c(1:6))
df_metadata <- subset(df_ST1T, select = c(1:6))
# Create bray-curtis distance matrix
dist_bray <- vegdist(df_abundance, method = "bray")
# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)
## Run 0 stress 0.04059044
## Run 1 stress 0.04059042
## ... New best solution
## ... Procrustes: rmse 8.989373e-05 max resid 0.0001311121
## ... Similar to previous best
## Run 2 stress 0.0405904
## ... New best solution
## ... Procrustes: rmse 1.546332e-05 max resid 3.034007e-05
## ... Similar to previous best
## Run 3 stress 0.2274369
## Run 4 stress 0.0405904
## ... Procrustes: rmse 1.309338e-05 max resid 3.302249e-05
## ... Similar to previous best
## Run 5 stress 0.04059044
## ... Procrustes: rmse 8.680016e-05 max resid 0.0001219727
## ... Similar to previous best
## Run 6 stress 0.1139412
## Run 7 stress 0.3314815
## Run 8 stress 0.117235
## Run 9 stress 0.04059063
## ... Procrustes: rmse 0.0001848809 max resid 0.0002566614
## ... Similar to previous best
## Run 10 stress 0.117235
## Run 11 stress 0.04059041
## ... Procrustes: rmse 2.967561e-05 max resid 4.274677e-05
## ... Similar to previous best
## Run 12 stress 0.1594032
## Run 13 stress 0.1139412
## Run 14 stress 0.04059041
## ... Procrustes: rmse 3.577167e-05 max resid 7.251582e-05
## ... Similar to previous best
## Run 15 stress 0.1139412
## Run 16 stress 0.04059054
## ... Procrustes: rmse 0.0001563901 max resid 0.0002213148
## ... Similar to previous best
## Run 17 stress 0.117235
## Run 18 stress 0.175473
## Run 19 stress 0.04059043
## ... Procrustes: rmse 5.695116e-05 max resid 8.059366e-05
## ... Similar to previous best
## Run 20 stress 0.04059044
## ... Procrustes: rmse 6.588413e-05 max resid 9.0608e-05
## ... Similar to previous best
## *** Solution reached
NMDS$stress
## [1] 0.0405904
# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## GenotypeSex 1 0.16453 0.19065 1.8845 0.1252
## Residual 8 0.69844 0.80935
## Total 9 0.86297 1.00000
# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.006147 0.0061474 0.3078 999 0.504
## Residuals 8 0.159771 0.0199713
Proximal small intestine figure
# Figures
#library packages
pacman::p_load(
ggthemes,
ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )
# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)
seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
by = c('GenotypeSex','Genotype','Sex'), sort = FALSE)
composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
geom_segment(data = seg,
mapping = aes(xend = oNMDS1, yend = oNMDS2),
size = 0.25) +
geom_point(data = centroids, size = 4) +
geom_point()+ #geo,_point()
stat_ellipse(linetype="longdash",size=0.3,color="gray32")+
coord_fixed()+
theme_bw() +
labs(title="")+
theme(legend.position = "bottom",
legend.title= element_text(size=10, face="plain"),
legend.text = element_text(size=14),
legend.key = element_rect(fill = "transparent", colour = "black"),
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 9), size=16),
axis.title.y=element_text(margin = margin(r = 3), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_x_continuous(limits = c(-1.5,1.5),breaks = seq(-1.5,1.5,0.6))+
scale_y_continuous(limits = c(-2.0,2.0),breaks = seq(-2.0,2.0,0.8))+
scale_color_manual(values=c("#30436d","#654e3a"))+
scale_fill_manual(values=c("#30436d","#654e3a"))+
scale_shape_manual(values=c(24,21,24,21))+
annotate("text",x=-1.5,y=-2,
label="Stress = 0.041",hjust=0, size=4)+
annotate("text",x=-1.5,y=2,
label="Proximal small intestine",hjust=0, size=5, fontface = "bold")+
annotate("text",x=-1.5,y=1.6,
label=expression(paste("Male (WT vs KO): ",italic("P") > 0.05)),hjust=0, size=4)
# Print
composition_GenotypeSex_final

Distal small intestine
set.seed(123)
# subset dateframe
df_abundance <- subset(df_ST2T, select = -c(1:6))
df_metadata <- subset(df_ST2T, select = c(1:6))
# Create bray-curtis distance matrix
dist_bray <- vegdist(df_abundance, method = "bray")
# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)
## Run 0 stress 0.04040786
## Run 1 stress 0.04040786
## ... New best solution
## ... Procrustes: rmse 8.780053e-06 max resid 1.891585e-05
## ... Similar to previous best
## Run 2 stress 0.1231312
## Run 3 stress 0.04040786
## ... New best solution
## ... Procrustes: rmse 3.686835e-06 max resid 7.196533e-06
## ... Similar to previous best
## Run 4 stress 0.04040786
## ... Procrustes: rmse 7.393811e-06 max resid 1.550481e-05
## ... Similar to previous best
## Run 5 stress 0.04040786
## ... Procrustes: rmse 8.841731e-06 max resid 1.950637e-05
## ... Similar to previous best
## Run 6 stress 0.04040786
## ... New best solution
## ... Procrustes: rmse 2.12053e-06 max resid 4.289854e-06
## ... Similar to previous best
## Run 7 stress 0.1231312
## Run 8 stress 0.04040787
## ... Procrustes: rmse 3.053269e-05 max resid 6.377193e-05
## ... Similar to previous best
## Run 9 stress 0.0472167
## Run 10 stress 0.1252416
## Run 11 stress 0.04040786
## ... Procrustes: rmse 1.870057e-06 max resid 4.302461e-06
## ... Similar to previous best
## Run 12 stress 0.04040786
## ... Procrustes: rmse 2.486584e-06 max resid 4.60433e-06
## ... Similar to previous best
## Run 13 stress 0.04040786
## ... Procrustes: rmse 2.887797e-06 max resid 5.481577e-06
## ... Similar to previous best
## Run 14 stress 0.04040786
## ... Procrustes: rmse 9.520293e-06 max resid 1.977447e-05
## ... Similar to previous best
## Run 15 stress 0.04040786
## ... Procrustes: rmse 2.542568e-06 max resid 4.669764e-06
## ... Similar to previous best
## Run 16 stress 0.0472167
## Run 17 stress 0.04040786
## ... Procrustes: rmse 3.249232e-06 max resid 5.789046e-06
## ... Similar to previous best
## Run 18 stress 0.1231313
## Run 19 stress 0.04040787
## ... Procrustes: rmse 1.450243e-05 max resid 2.672569e-05
## ... Similar to previous best
## Run 20 stress 0.1231313
## *** Solution reached
NMDS$stress
## [1] 0.04040786
# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## GenotypeSex 1 0.58745 0.54257 9.489 0.0071 **
## Residual 8 0.49527 0.45743
## Total 9 1.08272 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.030943 0.0309426 4.2672 999 0.082 .
## Residuals 8 0.058010 0.0072512
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Distal small intestine figure
# Figures
#library packages
pacman::p_load(
ggthemes,
ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )
# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)
seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
by = c('GenotypeSex','Genotype','Sex'), sort = FALSE)
composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
geom_segment(data = seg,
mapping = aes(xend = oNMDS1, yend = oNMDS2),
size = 0.25) +
geom_point(data = centroids, size = 4) +
geom_point()+ #geo,_point()
stat_ellipse(linetype="longdash",size=0.3,color="gray32")+
coord_fixed()+
theme_bw() +
labs(title="")+
theme(legend.position = "bottom",
legend.title= element_text(size=10, face="plain"),
legend.text = element_text(size=14),
legend.key = element_rect(fill = "transparent", colour = "black"),
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 9), size=16),
axis.title.y=element_text(margin = margin(r = 3), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_x_continuous(limits = c(-1.5,1.5),breaks = seq(-1.5,1.5,0.6))+
scale_y_continuous(limits = c(-2,2),breaks = seq(-2,2,0.8))+
scale_color_manual(values=c("#30436d","#654e3a"))+
scale_fill_manual(values=c("#30436d","#654e3a"))+
scale_shape_manual(values=c(24,21,24,21))+
annotate("text",x=-1.5,y=-2,
label="Stress = 0.040",hjust=0, size=4)+
annotate("text",x=-1.5,y=2,
label="Distal small intestine",,hjust=0, size=5, fontface = "bold")+
annotate("text",x=-1.5,y=1.6,
label=expression(paste("Male (WT vs KO): ",italic("P") == 0.0071)),hjust=0, size=4)
# Print
composition_GenotypeSex_final

Large intestine
set.seed(123)
# subset dateframe
df_abundance <- subset(df_LT, select = -c(1:6))
df_metadata <- subset(df_LT, select = c(1:6))
# Create bray-curtis distance matrix
dist_bray <- vegdist(df_abundance, method = "bray")
# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)
## Run 0 stress 0.119291
## Run 1 stress 0.119291
## ... Procrustes: rmse 9.908197e-05 max resid 0.0001913734
## ... Similar to previous best
## Run 2 stress 0.119291
## ... New best solution
## ... Procrustes: rmse 3.670601e-05 max resid 5.712139e-05
## ... Similar to previous best
## Run 3 stress 0.2385138
## Run 4 stress 0.119291
## ... Procrustes: rmse 0.0001347919 max resid 0.0002617218
## ... Similar to previous best
## Run 5 stress 0.119291
## ... Procrustes: rmse 1.54859e-05 max resid 2.286731e-05
## ... Similar to previous best
## Run 6 stress 0.2418937
## Run 7 stress 0.1192911
## ... Procrustes: rmse 6.660418e-05 max resid 0.0001270189
## ... Similar to previous best
## Run 8 stress 0.2825862
## Run 9 stress 0.119291
## ... Procrustes: rmse 4.138662e-07 max resid 7.575965e-07
## ... Similar to previous best
## Run 10 stress 0.1962946
## Run 11 stress 0.1192915
## ... Procrustes: rmse 0.0001413468 max resid 0.0002531473
## ... Similar to previous best
## Run 12 stress 0.2444607
## Run 13 stress 0.2532164
## Run 14 stress 0.119291
## ... Procrustes: rmse 6.078828e-05 max resid 0.0001113746
## ... Similar to previous best
## Run 15 stress 0.119291
## ... Procrustes: rmse 7.885862e-05 max resid 0.0001437482
## ... Similar to previous best
## Run 16 stress 0.119291
## ... Procrustes: rmse 7.812125e-05 max resid 0.0001434385
## ... Similar to previous best
## Run 17 stress 0.119291
## ... Procrustes: rmse 6.768069e-06 max resid 1.245338e-05
## ... Similar to previous best
## Run 18 stress 0.119291
## ... Procrustes: rmse 3.591897e-05 max resid 5.235938e-05
## ... Similar to previous best
## Run 19 stress 0.2518441
## Run 20 stress 0.119291
## ... New best solution
## ... Procrustes: rmse 1.59648e-05 max resid 3.229225e-05
## ... Similar to previous best
## *** Solution reached
NMDS$stress
## [1] 0.119291
# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## GenotypeSex 1 0.12799 0.1625 1.5522 0.1793
## Residual 8 0.65962 0.8375
## Total 9 0.78761 1.0000
# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0000402 0.00004020 0.0421 999 0.921
## Residuals 8 0.0076406 0.00095507
Large intestine figure
# Figures
#library packages
pacman::p_load(
ggthemes,
ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )
# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)
seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
by = c('GenotypeSex','Genotype','Sex'), sort = FALSE)
composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
geom_segment(data = seg,
mapping = aes(xend = oNMDS1, yend = oNMDS2),
size = 0.25) +
geom_point(data = centroids, size = 4) +
geom_point()+ #geo,_point()
stat_ellipse(linetype="longdash",size=0.3,color="gray32")+
coord_fixed()+
theme_bw() +
labs(title="")+
theme(legend.position = "bottom",
legend.title= element_text(size=10, face="plain"),
legend.text = element_text(size=14),
legend.key = element_rect(fill = "transparent", colour = "black"),
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 9), size=16),
axis.title.y=element_text(margin = margin(r = 3), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_x_continuous(limits = c(-1.5,1.5),breaks = seq(-1.5,1.5,0.6))+
scale_y_continuous(limits = c(-2,2),breaks = seq(-2,2,0.8))+
scale_color_manual(values=c("#30436d","#654e3a"))+
scale_fill_manual(values=c("#30436d","#654e3a"))+
scale_shape_manual(values=c(24,21,24,21))+
annotate("text",x=-1.5,y=-2,
label="Stress = 0.12",hjust=0, size=4)+
annotate("text",x=-1.5,y=2,
label="Large intestine",,hjust=0, size=5, fontface = "bold")+
annotate("text",x=-1.5,y=1.6,
label=expression(paste("Male (WT vs KO): ",italic("P") > 0.05)),hjust=0, size=4)
# Print
composition_GenotypeSex_final

ANOISM: Variance between groups and within groups
- ANOSIM analysis to test whether the variance between groups is larger than within groups (if it is larger, then P values is less than 0.05)
- Anoism analysis (R value from -1 to 1, R > 0 indicates between-group difference > within group difference; R < 0 indicates within group > between group difference; P value indicates the statistical significance)
# ST1T
## Define group
df_abundance <- subset(df_ST1T, select = -c(1:6))
df_metadata <- subset(df_ST1T, select = c(1:6))
## Bray-curtis
dist_bray <- vegdist(df_abundance, method = "bray")
## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)
##
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.1
## Significance: 0.174
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.184 0.252 0.324 0.332
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 14.00 26.0 34.0 45 25
## WT Male 3 18.25 22.0 34.5 44 10
## KO Male 7 9.25 14.5 29.5 37 10
plot(anosim.GenotypeSex_bray)

# ST2T
## Define group
df_abundance <- subset(df_ST2T, select = -c(1:6))
df_metadata <- subset(df_ST2T, select = c(1:6))
## Bray-curtis
dist_bray <- vegdist(df_abundance, method = "bray")
## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)
##
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.812
## Significance: 0.007
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.236 0.288 0.408 0.460
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 16 25.00 31.0 39.0 45 25
## WT Male 1 4.25 6.5 9.5 14 10
## KO Male 3 11.50 18.0 23.5 37 10
plot(anosim.GenotypeSex_bray)

# LT
## Define group
df_abundance <- subset(df_LT, select = -c(1:6))
df_metadata <- subset(df_LT, select = c(1:6))
## Bray-curtis
dist_bray <- vegdist(df_abundance, method = "bray")
## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)
##
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex)
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.12
## Significance: 0.768
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.203 0.284 0.308 0.392
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 9.00 18 39.00 45 25
## WT Male 3 17.75 29 34.75 38 10
## KO Male 8 21.25 24 28.50 33 10
plot(anosim.GenotypeSex_bray)
